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ABSTRACT 

The isorotation contours of the solar convective zone (SCZ) show three dis- 
tinct morphologies, corresponding to two boundary layers (inner and outer), and 
the bulk of the interior. Previous work has shown that the thermal wind equation 
together with informal arguments on the nature of convection in a rotating fluid 
could be used to deduce the shape of the isorotation surfaces in the bulk of the 
SCZ with great fidelity, and that the tachocline contours could also be described 
by relatively simple phenomenology. In this paper, we show that the form of 
these surfaces can be understood more broadly as a mathematical consequence 
of the thermal wind equation and a narrow convective shell. The analysis does 
not yield the angular velocity function directly, an additional surface boundary 
condition is required. But much can already be deduced without constructing 
the entire rotation profile. The mathematics may be combined with dynamical 
arguments put forth in previous works to the mutual benefit of each. An impor- 
tant element of our approach is to regard the constant angular velocity surfaces 
as an independent coordinate variable for what is termed the "residual entropy," 
a quantity that plays a key role in the equation of thermal wind balance. The 
difference between the dynamics of the bulk of the SCZ and the tachocline is due 
to a different functional form of the residual entropy in each region. We develop 
a unified theory for the rotational behavior of both the SCZ and the tachocline, 
using the solutions for the characteristics of the thermal wind equation. These 
characteristics are identical to the isorotation contours in the bulk of the SCZ, 
but the two deviate in the tachocline. The outer layer may be treated, at least 
descriptively, by similar mathematical techniques, but this region probably does 
not obey thermal wind balance. 
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1. Introduction 

Helioseismology has revealed a global pattern of differential rotation in the solar convec- 
tive zone (SCZ) comprising three distinct, readily apparent, regions. In the bulk of the SCZ 
(region 1), the isorotation contours appear mainly as quasi-radial spokes in meridional planes 
(corresponding in three dimensions to cones) . The contours change relatively little and only 
very gradually near the poles and equator, but in the boundary between the SCZ and the 
uniformly rotating radiative interior (the "tachocline" ) , the isorotation contours are abruptly 
wrenched into quasi-spherical shells (region 2). Finally, in the outermost few percent of the 
Sun's exterior layers, the SCZ isorotation spokes gradually acquire lengthening tangential 
tails stretching from the parent spoke toward the equator (region 3). The challenge for the 
theorist is to account for this global organization as economically as possible, a task that 
motivates the current paper. 

Despite the intrinsic dynamical complexity of a rotating, convecting fluid there are 
some indications that the solar rotation profile may obey rather simple mathematical laws. 
In particular, much of the convective zone may be governed by a reduced form of the vorticity 
equation, known as the thermal wind equation, or TWE (Thompson et al. 2003). The TWE 
involves a balance between inertial angular velocity gradients and baroclinic driving, the 
latter arising from angular entropy gradients that may be present. 

The TWE is not, by itself, sufficient to solve for the rotation profile, since it is a 
single equation containing two unknowns, the angular velocity fl and entropy S. It has 
been suggested, however, that S and Q may be a priori functionally related, and that this 
relationship could be exploited to solve the TWE. For example, Balbus (2009) pointed out 
that a certain class of magnetized baroclinic modes were marginally stable if constant S and 
fl surfaces coincided. With S a function of fl only, the TWE can be formally, but quite 
usefully, solved. The solution does not yield the entire spatial structure of fl, rather it yields 
a parameterized form for the surfaces on which O is a constant (see equation [3], below). 
The fit to the data is strikingly good. 

In a subsequent analysis, Balbus et al. (2009; hereafter BBLW) argued that magnetic 
fields were not essential to establishing the S — Q functional relationship. Under conditions 
expected in the SCZ, the effect of differential rotation on thermally convecting cells will be 
to distort surfaces of constant entropy into surfaces of constant angular velocity fl. More 
precisely, one expects the surfaces of constant residual entropy — the portion of the entropy 
in excess of the destabilizing radial entropy profile needed to drive convection — to tend 
to coincide with surfaces of constant angular velocity. In the TWE, a functional relation 
between residual entropy S' and fl works just as well as an 5" — Q functional relationship, 
since it is only the angular derivatives of S that enter into the equation: the difference 
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between S and S' is a function only of radius r. But the distinction between residual and 
total entropy is of course otherwise important. In numerical simulations of the SCZ that 
most closely resemble the Sun, for example, surfaces of constant S are nearly spherical shells, 
whereas surfaces of constant S' indeed correspond much more closely to those of constant 
angular velocity (Miesch, Brun, & Toomre 2006). 

The mathematical implementation of the constraint of shared surfaces of constant resid- 
ual entropy and rotation may be written 

S{r,e) = f{n') + Sr{r) (1) 

where S{r, 0) is the entropy, r and 6 are the standard spherical radius and colatitude angle, 
Vt is the angular velocity, / is an arbitrary function of and is the spherical entropy 
profile needed to drive and maintain convection. Equation ([T]) is then equivalent to the 
statement that the residual entropy S' = S — Sr constant on a surface of constant VP' . 

On a strictly mathematical level, however, Sr need not be a particular entropy profile — 
in equation ([T]), it evidently could be any function of r at all. Within the theory of BBLW, 
5*^. itself is never directly used in constructing the rotation profile. Instead, it is /(^^^) that is 
important. The arbitrainess of endows equation ([1]) with a much greater generality than 
is obvious at first glance, a point that we will develop in §2 of this work. 

If the ansatz ([T]) is used in the governing dynamical equation (the azimuthal component 
of the vorticity equation), a self-contained quasi-linear partial differential equation emerges 
for VL^ whose characteristics are precisely the isorotation curves of region 1. These contours 
take the very simple mathematical form 

r'^sin^e = A- B/r (2) 

where A and B are constant along a given curv^; this basic structure is in excellent agree- 
ment with the helioseismology data (BBLW). If the isorotation curve passes through some 
colatitude angle 9q at some radius r = tq, then the explicit equation for the isorotation 
curves is 

=in=9=(^^)'p»„-/3(^-l)]. (3) 

The dimensionless parameter /3 is a number of order unity, which might be different for 
different contours. 

It may be possible to develop this formalism to explain the structure of the isorotation 
contours even beyond the bulk interior of the SCZ. Recently, Balbus & Latter (2010, here- 
after BL) argued that a tachocline structure resembling that of the Sun emerges from the 



^More generally, sin'^ 9 = A + B^{r), where $ is the gravitational potential (Balbus & Weiss 2010). 
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assumption of a simple P2 spherical harmonic vorticity source that is present inside of some 
transition radius tt, and is otherwise r- independent. More precisely, if V/'Dr represents 
the radial derivative taken along a 6{r) characteristic given by equation ([3]), and if C is a 
constant, then the equation 



reproduces the SCZ isorotation contours for C = 0, and is able to closely describe the 
tachocline isorotation contours for C < 0. The angular velocity remains continuous at the 
transition radius r = r^; its derivative along the characteristic is discontinuous. In the 
bulk of the SCZ, the TWE characteristics and the isorotation contours are indistinguishable. 
Within the tachocline, they are very different curves. 

This development is quite suggestive, hinting that an understanding of the Sun's in- 
ternal differential rotation may not require a detailed dynamical knowledge of convective 
turbulence. Instead, the powerful constraint of thermal wind balance may be used to de- 
termine the angular velocity profile f2(r, 9) that is created by turbulent transport. At the 
same time, important questions arise which motivate the current work. In this paper, there 
is perhaps some progress toward elucidating some answers. But much remains to be done, 
particularly for the subsurface outer layers. 

Why is the match between mathematical solution and observations so good in the bulk 
of the SCZ? There seems to be a curious mismatch between the ability of the mathematical 
sin^ 6 = A — B/r curves to reproduce the helioseismology data so strikingly, and the more 
qualitative physical arguments that motivate the precise form of the thermal wind equation 
that is used. The theory seems to work too well. We would not expect, for example, constant 
residual entropy and constant angular velocity surfaces to match exactly. We would expect 
there to be a tendency to follow one another, and the most solar-like numerical simulations 
do show such a tendency (Miesch et al. 2006). But assuming an exact match leads to precise 
agreement between calculation and helioseismology observations. Is there something more 
going on here? 

Does thermal wind balance necessarily break down in the tachocline? BL did, in fact, 
interpret their equation this way, but is this essential? Is it possible to retain thermal wind 
balance with a different functional relation between angular velocity and residual entropy 
that keeps the goodness of fit intact? 

What is the physical explanation for the angular structure of the tachocline isorotation 
contours? Both the helioseismic morphology and the rotational dynamics seem to be com- 
pletely dominated by quadrupolar-like forcing. Why? For that matter, how accurately do 
we know that the angular behavior of the vorticity stress? Is it precisely quadrupolar? 
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What is happening in the Sun's outer layer? Do turbulent stresses become important 
(Miesch & Hindman 2011)? Can we guess the form of a simple ^-dependent forcing as was 
done for the tachocline? The surface distortion of the isorotation contours does not display 
P2 symmetry, it is strictly monotonic from pole to equator. The simplest possibility is sin^ 6 
forcing. How well does this work, and what is its physical origin? 

An outline of the paper is as follows. In §2, we review the results of earlier work (BBLW, 
BL) in which it was shown that thermal wind balance and informal arguments were used 
to construct precise mathematical solutions for the isorotation contours in the bulk of the 
SCZ through the tachocline. Wc then show that the all-important relationship between the 
residual entropy and must be true in a well-defined and physically interesting asymptotic 
limit, independently of other dynamical complications. These "complications," however, are 
likely to broaden the range of validity of the aymptotic limit. In §3, we construct a global, 
thermal wind model for the Sun's differential rotation profile, and show that it compares 
favorably with the hclioscismology data. The outer layers are treated on an ad hoc basis, 
with a deeper investigation of the issues presented in an Appendix. The change in the 
form of the isorotation contours apparent upon entering the tachocline does not indicate a 
breakdown of global thermal wind balance, but instead a local functional dependence of the 
residual entropy on both angular velocity and colatitude angle 9. In §4, we present a physical 
discussion of our results, explore its limitations, and suggest directions for future research. 



Throughout this paper, we will work with both spherical (r, 9, (p) and cylindrical (i?, 0, z) 
coordinates. As noted in the Introduction, r and 9 are respectively the spherical radius and 
colatitude angle. The quantity is the azimuthal angle, while R and z are the cylindrical 
radius and vertical Cartesian coordinate respectively. 

The governing dynamical equation of our analysis is the component of the vorticity 
equation. For steady fiows this is 



2. 



Thermal wind balance in the solcir convection zone 



2.1. Pure rotation model: governing equation 





or, since by mass conservation V'(pv) = 0, 
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where v is the velocity, u the vorticity, P the pressure, and p is the density. Equation ([6]) 
imphcitly defines a vorticity fiux within the divergence operator on its left side. The vector 
is of unit norm, in the azimuthal direction. 

The left side of equation ([5]) contains two terms, the first corresponding to advection by 
the poloidal velocity components of the vortensity (divided by cylindrical radius), the second 
to vorticity distortion by differential rotation. For a purely azimuthal velocity field, only the 
second (vorticity-distorting) term is present. The leads to what is commonly referred to as 
the thermal wind equation (Thompson et al. 2003): 

/aff\ sine /aff\ i 

An important simplification of equation (JTj) is effected by replacing, with impunity. In p by 
— (1/7) InPp"''', where 7 is the adiabatic index. Then, in terms of the entropy variable 
a = InPp^'^ and gravitational field g = —{l/p){dP/dr), we have: 

/dQ'^\ tane _ g /da\ 

\ dr J Q r \ 06 J J. ^r"^ sin 9 cos 9 \d9 J ^ 

On the right side, we have dropped the term proportional to {da/dr){dP/d9) in favor of 
{da /d9){dP/dr). This is justified because, while the radial gradient of a may well exceed 
its 9 gradient in the SCZ (region 1), because of entropy mixing by convection it does so by 
a factor far smaller than the pressure radial gradient exceeds its 9 gradient. We will adopt 
equation (|8]) as our governing equation, though it must break down, of course, inside the 
radiative zone (cf. §3). 



2.2. Establishing a functional relation between a and fi^ 

As discussed in the Introduction and in BBLW, the entropy variable cr appears in equa- 
tion dH]) only in the form of its 9 gradient. This means that any function differing from a 
by another function depending only upon r would serve just as well as a itself. There is, in 
other words, a gauge invariance in the proble 

This observation motivated BBLW to introduce the residual entropy a', defined by 

a'{r,9) = a{r,9)-ar{r) (9) 



^ A similar issue formally afflicts the angular velocity, which is insensitive to an additive function of 
R. This "geostrophic degeneracy" is lifted by fitting our solution to an explicit set of isrotation surfaces 
provided by observations. 
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where cr^ is the minimal adverse radial entropy profile needed to both drive and maintain 
ongoing convection. In this manner, we choose our gauge. The idea is that mixing of entropy 
would eliminate residual entropy gradients within a convective cell, but the resulting nearly 
constant value of a' within one convective cell need not be the same constant in another cell. 
BBLW went on to argue that the presence of shear would favor the tendency for long-lived 
coherent convective structures to lie also within constant Q surfaces. Under these conditions, 
constant a' and constant Q surfaces would tend to coincide. 

The mathematical implementation of this is a powerful analytic constraint. If a' = 
/(fi^), where / is a function to be determined, then 



(cf. equation [T]). In this view, the combination of convection and rotation in the SCZ 
constrains the two-dimensional entropy to be an additive function and r. 

Let us take a somewhat more formal viewpoint for the moment. On purely mathematical 
grounds, for this development to be valid, all that is required is that the constraint embodied 
by equation f lTU]) be satisfied for whatever reason. If it is, the ensuing thermal wind partial 
differential equation becomes 



where /'(f^^) is df /dQ"^, and we have made use of the P notation of equation (|4]). The solution 
to equation ( ITT]) is that Q is constant along the contours given explicitly by equation ([3]). 
It is remarkable that no knowledge of / is necessary to extract the basic functional form of 
the contours: is just a linear function of 1/r. Since these simple contours lie essentially 
on top of the helioseismology data (BBLW), the mathematical formalism at the very least 
constrains any deeper physical theory for the origin of Q{r,6). Equation (ITU]) and thermal 
wind balance are empirically true. The ultimate underlying physical cause remains to be 
settled, however: the physical argument laid out in BBLW may or may not be the correct 
one. 



In addition to a likely dynamical linkage between a' and Q, there is another simple 
feature of our problem that supports equation (ITUI) : the SCZ is rather thin. 

Imagine a Sun-like star with a slender outer convective zone. In the case of the SCZ 
itself, the central radius is about Tc = O.87r0 and our analysis extends iO.lr©, so the 



a{r,9) = f{n') + arir) 



(10) 




(11) 



2.3. Another approach to the a' — relationship 
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relative spread in the radial domain !S.rjrc is some 11.5%. The residual entropy o' is a two- 
dimensional function of r and 6', but it could equally well be locally regarded, by the implicit 
function theorem, as a function of r and ^ (the sign of should not matter j^l. We shall 
assume this is valid throughout the SCZ. Our motivation for choosing these two quantities 
as independent variables, as opposed to, e.g. angular momentum and r, is clear: we wish 
to turn the thermal wind equation into a simple mathematical constraint for the surfaces 
of constant f2. Had ^ not been one of our dependent variables, the task of extracting the 
form of the constant isorotation contours from the thermal wind equation would be much 
more difficult. Much the same technique is familiar in thermodynamics, where the choice 
of functions for independent variables is somewhat arbitrary, but often dictated by what 
quantity is being held constant during the transformation of interest. 

If the dispersion in ^ and r in the SCZ is such that o"'(fi^,r) never leaves the bilinear 
regime, then equation (fTOl) will be directly satisfied, with gauge invariance allowing us to use 
either a or a' . But an expansion in is not essential. If we expand cr'(fi^,r) in a Taylor 
series about r = rc, there obtains 



The second term, which we shall regard as a small, is reduced relative to the first by an 
explicit factor of 1 — r/rc ~ 0.1. But, as we argued in the Introduction, if in the dynamical 
regime of the SCZ there is already some tendency for constant a' and surfaces to coincide, 
which is another way of saying that any partial derivative of o' taken at constant Q will 
be small. The point is that "small" is by no means infinitesimal: the first-order correction 
to retaining the leading term in equation f|T2l) is quadratic, not linear in small quantities. 
Two 10% effects would combine to yield a 1% effect. It is the mutually reinforcement of 
the tendency of a' and Q surfaces to blend (even if it is not more than a trend), plus the 
narrow shell approximation |r — rd <^ rc, that renders the a'{r,6) = cr'(f2^) ansatz more 
robust than might at first sight seem apparent. The narrowness of the shell means that 
residual entropy variations are more beholden to changes in Q than to the spread in r, and 
the process of convection reinforces this trend by constraining the radial residual entropy 
gradients (at constant Q). The Sun happens to be in a felicitous theoretical regime from 
this point of view: it has both an order unity convective Rossby number (Miesch & Toomre 
2009) and a relatively narrow convective layer. 



■^There are of course some technical restrictions that accompany this remark. The variable fl^ should be 
well-behaved, free of internal extrema or saddle points over the domain of interest; surfaces of constant fl 
should not be spherical. These are not concerns for our present application. 
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Is it at all useful to consider the dependence a'{Q'^, 9)1 This seems counterproductive in 
the SCZ, where VL already has a dominant 9 dependence, and where we have just advanced 
arguments to the effect that a' should be considered a function of VL^ alone. But we can also 
turn this argument around. In the tachocline, VL is clearly dominated by its r dependence. 
Therefore, from a global viewpoint, adopting a a\Vt^,9) dependence would allow both the 
tachocline and the SCZ to be treated on the same footing, without necessarily abandoning 
thermal wind balance. In this formalism, a' goes from a two-dimensional function in the 
tachocline {VL^,9) to a simpler one- dimensional function in the SCZ (fi^). Such a dynamical 
structure would not be unknown in the study of rotating fluids: a classical two-dimensional 
Ekman boundary layer joining onto a one-dimensional interior Taylor column exhibits the 
same feature. In §3, this idea will be more fully developed. 



2.4. Intuitively understanding the convection zone 

2.4- i- Mathematical description 

The current approach has helped us to understand the answer to the first question 
posed in the Introduction: why do informal arguments work so well when our mathematical 
solutions are compared with observations? It is enlightening to examine the mathematical 
form of our simple SCZ solution with a view toward the helioseismology data. 

Let us start with the explicit solution for an isorotation contour in the form 

. 2n sin^^o + Z? /3 

sm 9 = ^ r (13) 

r y 

where y = t/tq. Then, expanding about some fiducial radius y, 



CI 



.2, /sin^^o+i? /3\ , .(W 2sin2^o + 2/3> 

^ = 3 -A^^V-V-)\-J + - (14) 



V yl VcJ \yt Vc 

which implies the existence of a constant 9 contour (a radial spoke) if 

1.5/3 
sin^ 9o + (3 

lies within the convection zone. The data show that in fact it does so for — 30°. For 
other angles not too near the poles or equator, the deviation from this class of contour is 
evidently not large. But neither is the alignment perfect: at larger 9q the spoke is canted 
more parallel to the z axis, while at smaller the cant is more orthogonal to the vertical 
direction. All this can be seen both in the helioseismology data, and the leading order term 
of a Taylor expansion of our simple solution. 



2.4-2. Dynamics 



What is special about Q ~ Q{6) that much of the convection zone is to first order well- 
described by this class of contour? It is by no means a direct consequence of our a' = 
ansatz, yet it is the most striking feature of the SCZ, and must involve the dynamics of 
heat convection and rotation. We defer an extended discussion of this important point 
to a subsequent paper (Schaan & Balbus 2011), but the basic picture is easily grasped. 
The coupled dynamics of entropy and angular velocity gradients are indeed at the heart of 
matter. The most efficient route for heat transport is spherical convection, radial motions 
being the most direct path to the surface. However, an inevitable secondary consequence 
of turbulent convection in an initially uniformly rotating system is the production of dif- 
ferential rotation, a sort of dynamical pollutant in this case. The presence of differential 
rotation introduces its own inertial force, pushing fluid parcels away from spherically radial 
motions toward cylindically radial motions. For a displacement ^, this acceleration is of the 
form — ^'VJl^. (Note that this is not the Coriohs force.) If spherical radial displacements 
dominate, this oblateness- inducing force is formally eliminated when Q = Q.{6). Radially 
convecting elements could then not help but move in constant Vt surfaces. This closes the 
circle of theoretical self-consistency: the gross pattern of differential rotation in the Sun 
looks the way it does because in the presence of unavoidable build-up of differential rotation, 
VL ~ Vt{9) would produce a pattern that would not interfere with the most efficient form 
of heat transport. The top priority of the convective zone would remain uncompromised. 
Moreover, we have in the previous section noted that the description of the curves is 
not precisely ^2(6*); the curves are canted slightly poleward everywhere but at the highest 
latitudes. This slight poleward cant emerges, in fact, from the linear theory of convective 
displacements (Schaan & Balbus 2011) as well, and ^'VJl^ remains small. 



3. Global solutions of the isorotation contours 

3.1. Entering the tachocline 

As one crosses from the SCZ into the tachochne, the isorotation contours change their 
character abruptly. It is possible that this is due to the increased importance of additional 
dynamics beyond pure rotational motion, and the thermal wind equation breaks down. On 
the other hand, there is little reason to assume that the precise a' — fi^ relation in the SCZ 
should hold here; the r component of Vr2 becomes vastly more important in the tachocline. 
Let us assume that thermal wind balance continues to hold in the tachocline. We will 
make use of the ansatz a' — a'{fl'^,9), which is particularly apt, since ft now carries a 
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strong r-dependence. This pronounced dependence on radius is of course a consequence of 
the boundary condition enforcing uniform rotation at the radiative zone interface. Once 
this zone is breached, the radial entropy gradient starts to rise, and the approximation of 
neglecting the term proportional to ds/dr in the baroclinic contribution to the thermal wind 
equation eventually breaks down. But by this point, Q will be so close to uniform rotation 
that the shift of the isorotation contours is unlikely to be important. 

We turn now to an investigation of thermal wind balance in the tachocline. 

3.2. The ansatz a' = 9) 

The thermal wind equation acquires a very different mathematical structure if a' is taken 
to be a function of and 6, as opposed to and r. We will assume that this functional 
dependence is well-defined locally in r for the tachocline region. The, the equation formally 
becomes 

on'' 

dr 

This should be contrasted with the case a' = a'(fi^,r), for which the equation takes the 
form 

dr 

What is striking here is that equation (fT6|) has precisely the same mathematical form invoked 
by BL to account for the rotation pattern in the tachocline. (The computed rotation pattern 
fit the data well.) But this earlier study viewed the right side of the equation as arising 
from externa/ vorticity stresses. Here, we see that this same equation could equally well arise 
without the need to appeal to external dynamics, emerging instead within the context of 
thermal wind balance. 

Notice the important distinction between da'/dQ'^ at constant 6 versus constant r. At 
constant 6, the derivative is perfectly well behaved in the tachocline. At constant r, the 
derivative is ill-behaved: a significant change in a' can be accompanied by virtually no 
change in Q. For this reason, equation ( !T6|) is to be preferred. 

Comparison of equations (1161) and (I17p reveals another simple but important fact. If a' 
is a function of fl^ alone in the SCZ (and it certainly appears to be so), equation (ITB]) is valid 
globally. The partial 6 derivative on the right vanishes in the SCZ, and comes alive within 
the tachocline. At the very least, this is a useful organizational framework for understanding 
solar rotation. The truly daunting problem of mastering the underlying dynamics of the 
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tachocline can be decoupled from the strictly mathematical problem of finding self-consistent 
solutions to equation f|T6|) . For the latter, we may appeal directly to the observations, and to 
very general and simple dynamical constraints. While not directly addressing the underlying 
physical structure of the tachocline, this approach is not without interest: a compelling fit 
would be strong evidence for the assumption of thermal wind balance. 



3.3. Tachocline solution 

In BL, the characteristics of the SCZ continued into the tachocline, and the entire inho- 
mogeneous term on the right side of the governing equation was assumed to be proportional 
to a function of 9 only, the spherical harmonic P2{cos9). 

In the current theoretical formalism, we have created somewhat tighter mathematical 
constraints for ourselves. The SCZ characteristics can now extend into the tachocline (which 
BL found produced good agreement with the data) only if {da' /dfl'^)e is a function of 
alone. This means that {da' / 89)^2 can be a function only of 9. This in turn means that the 
right side inhomogeneous term of equation f|T6|) has a steep g/r"^ ~ l/r"^ radial dependence. 
This radial dependence was not present in BL, since the entire transition zone (from O.77r0 
to O.65r0 was viewed as a local boundary layer. How does the additional radial dependence, 
a new feature of the current model, affect the shape of the tachocline isorotation contours? 

In the tachocline, equation (fT6|) may be written 

_ 9 (da'\ _ 

Vr sin 9 cos 9 \d9 J ^2 

The 9 dependence of the right side of the equation may be inferred from helioseismology 
data as being approximately proportional to —P2{cos9), or sin^ 9 — 2/3. We will consider a 
slightly more general form, sin^ 9 — e, allowing e to be determined by the best fit. (Note that 
the 9 dependence of a' itself must then involve a superposition of P4(cos^) and P2{cos9), as 
found in the simulations of Miesch et al. [2006].) Consider the equation 

^ = -,{sin^9-e). (19) 

The formal solution to equation (fT9|) is that VL^{r) is given along a characteristic curve by 

l]2(r) = fi2(ri) + I" ^[sin^ 9 {r)-e]dr (20) 

Jri ^ 



with 

sin^ I 



{^y[s,,e„-,(^}-i)]. (21) 



- 13 - 



The coordinates tq and represent the starting and defining point of the characteristic (the 
effective surface), and ri is a fiducial radius marking here the starting point of the tachochne. 
In general, we will use equation fl2U]) with r < ri, reversing the integration limits and changing 
the sign. The /3 parameter in general will depend upon 6*0, though it is simplest to take (3 
to be a global constant, and we will often do so. The radius ro will always lie outside of the 
tachocline on the constant-il portion of the characteristic. We use the notation ri = r-p for 
start of the SCZ-tachocline boundary, and we may take fi^(ri) = f2Q(sin^ 9o). Then, along a 
characteristic (12 ip . the variation of angular velocity is given by 

n^ir) = fi2(sin2 Oo) + eCF3 - Crl{P + sin^ ^0)^5 + C/^r^Fg (22) 

where 



To apply this with a minimum of mathematical overhead to an interesting case, consider 
a linear form for nQ(sin^ ^o), 

fi2(sin2 ^o) = nl + Ql sin^ ^0, (24) 
with and fi^ constants. This may be regarded as a simplified model of the Sun's surface. 

For this particular form of fio(sin^^o); the solution for Q'^{r,9) is found by eliminating 
sin^ 9o between equations and (12^ . The equation for the isorotational contours is then 
found to be 



A - ecir^Fs + (3cirlF^ - /^cirgFg ^ ^ ^ _ ^ 



1 - cirgFs 

where we define the dimensionless constants A and ci by 



(25) 



^ = ^F^' ^^ = 7?F ^^^^ 

Note that the contour-labeling collection of constants comprising the quantity A is numer- 
ically the same as sin^ 6*0. By choosing ci and e appropriately, equation ( 125|) represents a 
global solution for both the bulk of the SCZ and the tachocline. 



3.4. The outer layer 



In the outermost layer of the Sun, r > 0.96rQ, the convection is likely to be vigorous, 
with velocities comparable to the sound speed. In this zone, there is no reason to expect 
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that a simple thermal wind balance prevails in the vorticity equation, and good reason to 
suspect that turbulent forcing is involved (Miesch & Hindman 2011). For deriving the form 
of the isorotation contours however, it is useful to have a simple phenomenological descrip- 
tion of this region. An inspection of the helioseismology data indicates a sin^ 9 rotational 
morphology in the subsurface layers. We have accordingly found that setting 

mf_ ^ Co sin' e 

Vr Tq — r 

works extremely well as a governing equation for the outer layer isorotation contours. The 
right side of this equation is a departure from thermal wind balance, and may be the outcome 
of some combination of convectively-induced Reynolds stresses, meridional circulation, and 
magnetic braking. Here, as earlier, V/Vr is the radial derivative along the SCZ contour, 
and Co is a universal constant. In the formalism of Miesch & Hindman (2011) the baroclinic 
term is dropped, whereas here it is incorporated in the V/Vr operator. The right side 
of equation (!27j) would then correspond to the Q term of Miesch & Hindman, divided by 
RcosO. Whatever the origin of the right side forcing of (I27p . this equation may be solved 
for the constant (3 case by characteristic techniques entirely analogous to those used for the 
tachocline solution f l25|) . In fact, the equations are simpler, since r may be set equal to 
everywhere, save the expression r — Tq. 

An investigation of the form of equation ( 1271) based on retaining poloidal circulation 
terms in the vorticity conservation equation (|5]) is presented in the Appendix. We find 
that the sin^ 9 dependence can be readily reproduced in such a picture, but the singular 
l/(r0 — r) behavior seems to require additional dynamics. It is possible that the increasingly 
steep radial entropy gradient plays a role via the as yet neglected baroclinic term proportional 
to the product of dP/dO and da /dr. 



3.5. Results 

3.5.1. Constant and varying (5 

The analytic solutions that extend our characteristic techniques into the tachocline 
and SCZ boundaries are rigorous, we reemphasize, only if the /3 parameter is a universal 
constant. Otherwise there is a nontrivial coupling between the characteristic equations for 
dO/dr and d^ jdr in the boundary layers. (But not, it must be stressed, in the bulk of the 
SCZ, which may handled rigorously in the case of a varying /3.) We will develop accurate 
numerical models for this case in a later study, but for the present we will begin with 
constant /3 models. These work well enough (remarkably so given their simplicity) to foster 
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confidence in our fundamental approacli. We will then take our expressions for the constant 
/9 solutions, and in effect abuse them by using the same expressions with a varying /3 chosen 
to match the contours in the bulk of the SCZ. Ignoring the differential equation couplings this 
would introduce with the tachocline is justified if such couplings were small, here a marginal 
assumption. But the effect of this simple expediency on the boundary layer isorotation curves 
produces an extremely striking match to the data (cf. Figure [5]). This strongly motivates 
pursuing a more rigorous numerical solution of the governing partial differential equation, a 
project we defer to a later publication. 

Figures (1-3) show the constant (3 isorotation contours (equation [25]) for (3 = 0.55, 
0.84, and 1.1 respectively. The first value is chosen to match high latitude structure, the 
second to match the diverging tachocline contours at ^ 60°, and the final /3 value to match 
the equatorial contour structure. Overall, the global constant /3 models clearly resemble the 
helioseismology data of figure 4. In particular, the large /3 contours near the bifurcation 
point of the tachocline display a "plateau" followed by a precipitate, cliff-like structure that 
is evident in the data. 

In figure (5) we have overlayed the data with an analytic solution. To effect this, we 
have used equation fl25|) but with /3 a function of sin 6*0: 

/3 = 2.5 sin^ ^0-2. 113 sin 00 + .8205. (28) 

The overlay shows an excellent, detailed, global match. 



3.5.2. What IS e? 

In figures (1-3), e has been set to 0.75. This parameter sets the location at which the 
isorotation contours bifurcate in the tachocline. What sets its value? 

BL set e equal to 2/3 on the basis that this fit the data rather well, and because 
sin^ 9 — 2/3 was proportional to the spherical harmonic P2[cos{9)], a natural leading order 
structural modification to spherical symmetry. But a closer investigation shows that there 
is a better data fit. In figures (5) and (6) we show overlay plots with e = 2/3 and 4/5 
respectively. While these fits are reasonable, they are clearly inferior to e = 3/4. Is there 
something special about this value of e? 

One possible explanation is dynamical. For example, what is the torque exerted by 
our solution on the inner, uniformly rotating radiative zone? We assume that the couple is 
proportional to Rdfl/dr, the r(j) component of a viscous stress tensor. We denote 

n + 1 , , 

en = —— 29 

n + 2 
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so that ei, 62, and 63 correspond to our three choices. Then the following elementary integral 
identity is satisfied: 



The case n = 3 is relevant to the calculation of the torque exerted by the tachocline on the 
inner uniformly rotating core since this torque is directly proportional to the integral on the 
left side of the equation. (One factor of sin 6 arises from the r0 stress tensor component, 
another from the lever arm rsin^, and the third from the spherical area element.) In other 
words, when e — — 0.8, the torque exerted by the tachocline on the interior vanishes. 
Since the couple to the core must act on time scales no shorter than evolutionary, we might 
therefore expect e = 0.8 to yield the best fit. 

The best fit, however, is clearly e — e2 — 0.75. This is close enough to 0.8 that there may 
be something correct about our argument, but far enough off the mark that there is surely a 
piece of the puzzle missing. Moreover, it really is the case that e = 3/4 is the best fit in our 
models, not just among the three values we display here. While there is no compelling reason 
to believe that the value of e may be accurately deduced from the idealized assumption of 
a vanishing viscous torque, it is puzzling that an equally simple, but quite distinguishable, 
integral constraint seems to fit the data so well. It is as though we had mistakenly included 
an additional factor of sin 6 inside our viscous torque integral. 

For the present, the explanation for why e = €2 = 0.75 and not €3 = 0.8, corresponding 
to a true vanishing stress at colatitude 6 = 60° rather than 3.4° further south, remains 
elusive. 



In this paper, we have investigated global models of the angular rotation profile of the 
solar convective zone (SCZ) using the toroidal component of the vorticity equation and under 
the assumption of a dominant thermal wind balance. The model is global in that it seeks 
to unite both the tachocline and the bulk of the SCZ within the thermal wind formalism. 
The subsurface outer layers have also been treated at a quantitative, but more descriptive, 
level. In effect, we have modeled the subsurface outer layers by an ad hoc source term to 
the thermal wind equation. 

The general nature of the isorotation surfaces in the joint SCZ-tachochne structure 
seems to be a mathematical consequence of a few simple principles, each noncontroversial 
by itself: i) thermal wind balance, as noted above; ii) a narrow SCZ which may be treated 
locally in r; iii) a tendency for surfaces of constant angular velocity and residual entropy (see 




(30) 



4. Summciry 
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below) to blend; iv) an inner boundary condition of uniform rotation at the tachocline base; 
v) vigorous radial transport in the bulk of the SCZ; vi) a small couple between the tachocline 
and the core. When modeled by an idealized viscous prescription and the simplest possible 
deiviation from spherical symmetry, assumption vi gives close, but not perfect, agreement 
with helioseismology data. 

As in earlier studies (e.g. BBLW, BL), we have introduced and exploited the concept 
of residual entropy, in essence the true entropy minus a suitable average over angles. The 
residual entropy is, in fact, defined only up to an additive function of r. Numerical simula- 
tions are particularly well suited to this concept: a Schwarzschild-unstable, user-determined 
entropy profile S{r) is imposed to drive the requisite convective turbulence, and what one 
may call the "response" entropy S{r, 6) is then the quantity that is calculated (e.g. Miesch & 
Toomre 2009 and references therein). While S may formally differ from the residual entropy 
S' by compensating the radial background S{r), in practice it is not diffiult to extract a 
suitable S' from the numerical simulations. 

The heart of our approach is to regard the two-dimensional residual entropy as a function 
of Vt^ and 9 in the tachocline, but a function of and r in the bulk of the SCZ. In §2.3 
we have argued that the r dependence of a' within the SCZ is ultimately unimportant. In 
fact, it seems possible to arrive at the appropriate form of the SCZ thermal wind equation 
for this latter case via multiple routes; the advantages of theoretical models that lead to 
the coincidence of constant a' and VL surfaces have been discussed elsewhere (Balbus 2009, 
BBLW, Schaan & Balbus 2011). In any case, such behavior enjoys very clear numerical 
support (Miesch et al. 2006). 

The explicit bulk-SCZ and tachocline solutions, equations (l2Ti) and (!25|) respectively, 
work reasonably well even in the constant /3 limit, and strikingly well when /3 is permitted 
to vary. Only the details of the mathematical solutions, not their essential structure, are 
sensitive to the choice of free parameters. At this point there seems little doubt that thermal 
wind balance is, in fact, the overall organizing principle for the tachocline and interior SCZ. 

We are, however, far from a complete solution. The mathematical equations contain free 
parameters, which must be chosen to fit the helioseismology data. We lack an understanding 
of how they are determined from first principles. For example, the existence of a tachocline 
"bifurcation point," where the isorotation contours diverge either toward the pole or equator, 
is clearly a consequence of the ^-independent boundary condition imposed by the uniformly 
rotating core. The location of this bifurcation point depends upon our e parameter, which 
reproduces the data well if e = 0.75, and much less well if e changes even by a little. Why 
e should be close to, but distinct from, the value of 0.8 (corresponding to vanishing viscous 
torque on the core) is not yet clear. Finally, the outer layers are not well understood. 
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Once again, the mathematical description of the isorotation contours seems simpler than 
the underlying physics that gives rise to them. A good match to the helioseismology data 
is obtained from an approach that uses thermal wind balance plus a simple sin^^/(r — r©) 
external source. 

We may thus end on a positive note. If the work presented here is correct in its essen- 
tials, then the beginnings of a deeper dynamical understanding of the rotation of the solar 
convective zone are at hand. It is gratifying that there is a certain mathematical inevitabil- 
ity of the gross features of the Sun's internal rotation pattern following from thermal wind 
balance and a tight coupling between the rotation and the residual entropy isosurfaces. It 
is likewise encouraging that the more complex appearance of the isorotation contours in the 
tachocline and outer layers seems to emerge from relatively simple mathematics. In short, 
there is no reason to think that the dynamics of the solar rotation problem is intractable. 
The blend of numerical simulation and analysis promises to be a powerful combination for 
elucidating the rotation profile of the convective regions of the Sun and of other types of 
stars. 
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Appendix: the outer layers 

The rotation pattern of the outer layers is characterized by a robust sin^ 6 angular 
pattern, and the stress associated with this seems to depend on radius as l/(r — Rq). We 
assume that there are large scale poloidal velocity components in the outer layers, arising 
perhaps from turbulent stresses present in this region (Miesch & Hindman 2011). The 
introduction of nonuniform poloidal velocities Vr and ve into the flow causes the (f) component 
of the vorticity, a;^, to come into existence. Because of our assumption of axisymmetry 
{d/d(f) — 0) there are no changes at all to the poloidal components of uj. The addition of Vr 
and vg velocities also do not cause any cross terms between Q and the poloidal components in 
the azimuthal vorticity equation. This means that the poloidal components of the velocity, 
which advect changes in a;^, behave as "external forcing" with regard to the angular velocity 
fl in the governing partial differential equation. 

Start with the simplest possible circulation pattern consistent with quadrupolar sym- 
metry, 

V0 — A{r) sin 9 cos 9. (31) 

Since mass conservation implies V'(pv) = 0, the angular dependence of Vr must be the same 
as 

= — (2cos — sm fc/j = — (3cos e/ — 1) (32) 



rsin^ 89 r r 

In other words, Vr has the angular dependence of P2(cos^). Letting 

Vr ^ B{r) {3 cos'^ 9 -1) (33) 

we find that B and A are related to one another by the mass conservation equation 

I^AEI^ = -rA (34) 
p dr 



This preprint was prepared with the A AS IM]eX macros v5.2. 
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Next, we evaluate the terms in pv •'V {u ^ / pr sin 6) from the vorticity equation ([5]). The 
terms involving the operator Vrd/dr are 



pvr 



d 



d{rve) 



dvr 



dr ipr^ sin 9 dr pr"^ sin 9 89 
which simplifies to 



pVr COS 9 



d 

dr 



1 d{rA) QB 



pB{3cos^9- l)cose 



d_ 

dr 



1 d{rA) 6B 
pr"^ dr pr'^ 



(35) 



(36) 



This leads to a forcing term that is directly proportional to P2(cos^), which is not what is 
observed. By contrast, the terms involving V0d/d9 are 



— —cos 9 
r^ 09 



d{rA) 
dr 



+ 6B 



A sin^ 9 cos 9 



d{rA) 
dr 



+ QB 



(37) 



This gives a second forcing term proportional to sin^ 9, which is consistent with the obser- 
vations. 



Gathering the two terms together, 



pv-'V 



UJ4, 



pcos9 



pr sin 9 

which, on the right side, is equal to 



B{3cos'9-l)- 



d Asin^6' 



dr 



1 d{rA) 6B 



pr"^ dr 



+ 



pr^ 



(3^ 



p cos 9 



r^/o 2/1 . \ ■ 2n ^ dipr^B) 

Bi?,cos^9-l)— + sm^9- ' 
dr 



pr"^ dr 



1 did, 9 ^, QB 
pr^ dr p dr pr"^ 



(39) 



Thus, there are two possible angular behaviors of the vorticity advection terms that 
emerge from the analysis: P2(cos^) and, possibly in superposition, sin^ ^. For the outer 
layers, however, we require forcing that is proportional only to sin^ 9. This means that the 
terms in the second square bracket of equation (139|) must be a constant, which we denote 
as C. Standard models for the outermost layers in the Sun seem to give something close to 
p = Pq{1 — xY ioT the density, where po is a characteristic density and x = t/Pq, the radius 
in units of the solar radius. Then, with u = 1 — x, p = pou^, this constancy condition may 
be written as an inhomogeneous differential equation for B: 



did 



du V? du 

The solution for the homogeneous equation is 



{v!^B) + 6B = CRIpou^ 



B [homo) = Ci h C2- 



u 



u 



(40) 



(41) 
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and the inhomogeneous solution is (to leading order in small u) 

B [inhomo) = Cpo^| (42) 

32 

The inhomogeneous solution is well behaved near the solar surface (i.e., for small u), 
and the I2 solution is well-behaved at m = 0, but the K2 solution is singular. The mass 
flux, however, is finite. Moreover, the disturbance should disappear as we move into the 
convective zone, so that our solution must decrease with growing u. This is not compatible 
with I2. We therefore retain the K2 solution — indeed, we are looking for singular velocity 
behavior near the surface. Then, for small u we have 

_ 1 djpB) _ 1 dju^B) ^^g^ 
p dx du 

or, 

c^ d{u'K2{y/Qu)) _ Cl d{U^K2{U)) 

du u^U dU ^ ' 

where U = a/6m. Using a standard Bessel function identity, this becomes 

A = -^U'K.iU) = -^K,iV6u) = -^K^iVGu) (45) 

Therefore, to extract a solution with sin^ 6 angular forcing in the outer layers from the 
assumed circulation pattern, the quantity uj^/ {pr saiO) and radial mass flux must both be 
constant with depth, and the 6 velocity sharply decreasing with depth. 

If we now pull together the results of this mathematical excursion and return to equation 
(l38l) to determine the nature of the forcing we find i) it is rather straightforward to produce 
sin^ 9 forcing; and ii) the sign of the forcing can be chosen by a suitable choice of C and 
Cl. These are perhaps two important successes. But in spite of the singular behavior of the 
K2 function near the surface, in this model the vorticity stress is not singular. The regular 
behavior is due to the appearance of the surface density p in the equation, which vanishes 
at the surface. This is an important shortcoming of the model, as the data seem to require 
logarithmic "quasi-singular" behavior of the constant Vt contours near the surface. It may be 
that this feature can be captured by a more careful consideration of the radially steepening 
entropy gradient, which has been neglected here, lying outside the scope of our analysis. 
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Q. Contours 




X / r 

sun 

Fig. 1. — The isorotation contours of the Sun according to equation ( 125|) with (3 = 0.55, 
ro = 0.896, n = 0.96, tt = 0.77 (solar radii), h = 1.745. 

D. Contours 




X / r 

sun 



Fig. 2.— As in Fig. [1] with /3 



= 0.84. 




Fig. 4. — Global Oscillation Network Group (GONG) isorotation contours (courtesy R. 
Howe). 
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Q. Contours 




Fig. 5. — GONG data (black) overlayed with fit from equations fj25|) and fl28|l (red), e = 0.75. 
Other parameters as in Fig. [1]. The contours near the point of bifurcation at the radiative 
core are very well-modeled. 



- 25 - 



Q. Contours 




Fig. 6. — GONG data (black) overlayed with fit from equations fl25|) and fl28|) (red), e = 2/3. 
Other parameters as in Fig. [1]. 
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Q. Contours 




Fig. 7. — GONG data (black) overlayed with fit from equations fj25|) and fl28|) (red), e = 0.8. 
Otlier parameters as in Fig. [1]. 



